Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Feature force lever #157

Merged
merged 47 commits into from
Feb 9, 2024
Merged

Feature force lever #157

merged 47 commits into from
Feb 9, 2024

Conversation

chryswoods
Copy link
Contributor

Changes proposed in this pull request:

This PR introduces two main pieces of code:

  1. Support for LambdaSchedule equations to be specialised for individual levers for individual forces, and, indeed, for individual perturbable molecules. This will let us create schedules for absolute binding free energy calculations. This is described in the new Part 7 of the tutorial, which is leading up to absolute binding free energies (we will add more when we bring in the latest Boresch restraint code in that PR, plus add in some simple functions for creating annihilating merge molecules etc.)

  2. Support for better constraints in OpenMM - this includes dynamic constraints, the "-not-heavy" versions of the other constraints, and lots of small bug fixes and tweaks so that we get much better agreement with somd1.

There's also lots of other little bits of functionality as described in the changelog.

  • I confirm that I have merged the latest version of devel into this branch before issuing this pull request (e.g. by running git pull origin devel): [y]
  • I confirm that I have added a test for any new functionality in this pull request: [y]
  • I confirm that I have added documentation (e.g. a new tutorial page or detailed guide) for any new functionality in this pull request: [y]
  • I confirm that I have added a changelog entry to the changelog (we will add a link to this PR as part of the review): [y/n]
  • I confirm that I have permission to release this code under the GPL3 license: [y]

Suggested reviewers

@lohedges, @mb2055

Any additional context of information?

There is a breaking API change for LambdaSchedule, as described in the changelog. The tutorial and tests have been updated to the new API.

…ng them to the

documentation and options objects. Just running some simulations now to make sure
they are working as expected. Debugging showed the right bonds and angles were
not being constrained.

I've updated the defaults so that ``h-bonds-not-perturbed`` etc are the
default constraints if none are set. I've also added a discussion of PME
versus RF to the tutorial, and updated the example FEP script to
use ``h-bonds-not-perturbed`` and ``RF`` by default.

On my computer, the entire methane/ethanol FEP simulation now takes
about 3 minutes :-)
adding tests to check that the constraints are correct

WIP
different schedules for different molecules (eventually will be parts
of perturbable molecules)

This will let us use a single lambda value to control perturbations
of multiple molecules (or parts of molecules) that perturb their
parameters at different rates in different forces (or even go in
different directions).

This will let us do absolute binding calculations where we only turn
off the intermolecular energies (switch off ghost/non-ghost but keep
on ghost/ghost), or to do QM/MM end state corrections
(switch off ghost/ghost while keeping ghost/non-ghost the same)

I need to add tests and documentation
…l correct

if we link to the perturbed properties and/or swap the end states
…harge correction

in the ghost forces. This will be used to control that correction, e.g. so that
we can separate the ghost/ghost coulomb force from the ghost/non-ghost, e.g.
in absolute binding calculations where we don't want to decouple the ligand
from itself.
Fixed a bug where the h-bonds constraint would look at the highest mass
of perturbed atoms, when it needs to look at the lightest mass
…X), where

X can be the number of steps or the time between updates.

Also fixed tests_openmm_constraints.py
same energies. Also fixed a small bug where setting the dynamics
schedule didn't trigger a change in energy if lambda doesn't change
…dynamic_constraints"

option to mol.dynamics()
…ions to be

counted twice. They were added to the Custom force used for the perturbable torsions,
but then also added to the PeriodicTorsionForce for non-perturbable torsions.

This is because the incorrect mirrored form of the DihedralID was added to the
…aframe

with energies in reduced units. Also, the sr.morph.to_alchemlyb function
will preserve attributes for the combined dataframe
… hand-written ones

This was more painful than I expected... But all working now
…ile and the PerturbableOpenMMMolecule reporting

functions working. Sorting out the function call layout for a tutorial, and also adding unit tests

WIP
…n in the OpenMM context

You can get these via

```
p = mol.perturbation().to_openmm(constraint="...", etc)

print(p.get_changed_atoms())
print(p.get_changed_bonds())
print(p.get_changed_angles())
print(p.get_changed_torsions())
print(p.get_changed_exceptions())
print(p.get_changed_constraints())
```

I have also added a function to zero torsions in the end state where
any atom is a ghost. This is `sire.morph.zero_dummy_torsions`.
…use the

new functionality. Now working this up, updating the functionality as needed
to make a nice tutorial.

There are some changes to the code which I've described in the changelog updates.

In particular, the behaviour of ``mol.perturbation().link_to_reference()``
has changed, meaning that you don't need to add a ``.commit()`` on the
end (indeed, it will now give an error).

The right fix for this is to use the new, higher level functions to
link to the reference or perturbed states, e.g.

```
mols = sr.morph.link_to_reference(mols)
```
…onstraints "X",

but will not constrain perturbing bonds *unless* they involve a hydrogen in
any end state.

I have switched this to be the default constraint type, as it is the closest match
to somd1. The real match would be
`constraint="h-bonds", perturbable_constraint="bonds-not-heavy-perturbed"`

I have also improved the code that detects "light" atoms (i.e. hydrogens).

This now uses the maximum mass in any end state being below 2.5, or if
an element in any end state is a H, or if the ambertype in any end state
matches "H*" (case-insensitive).

These three ways of matching are all switched on by default. They can
be disabled via options described in the OpenMM cheatsheet.
…nd also

how to set it for specific perturbable molecules.

Also added a record in the schedule for the name of all of the forces,
so these are easier to query.
I've also exposed an option to control the MC barostat frequency
…softening potential.

Set the default to 1 A, so that this matches the somd1 value (and interpretation).

The previous value was equivalent to 10 A.
Also updated the somd coulomb constant to be more precise, matching that
used in sire/openmm. This doesn't change any energies, but I am happy
having a good level of consistency.
Also BREAKING CHANGE in the LambdaSchedule API, as described in the changelog
@chryswoods
Copy link
Contributor Author

I've cancelled the PR as fixing the Python 3.12 issue was quicker than I expected. That PR also includes this PR, so CI on that will also test this (trying to save GH Actions time).

@lohedges
Copy link
Contributor

lohedges commented Feb 9, 2024

Thanks for this. I've been using this branch locally for a few weeks now so don't foresee any issues. I'll re-run the BioSimSpace and somd2 tests locally to make sure all is okay. I think the biggest issue will be moving my QM/MM tutorial once I merge into the feature_emle branch ;-)

@chryswoods chryswoods merged commit 7329ccc into devel Feb 9, 2024
0 of 5 checks passed
@chryswoods chryswoods deleted the feature_force_lever branch February 9, 2024 18:13
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

[BUG?] Does the alanine dipeptide test system have a mixture of combining rules?
2 participants